function WayPoint3D()

global n WP3D COG npoints Resolution
global WP3DStandard TimeResolution
global pointIdentifiers f1 f2

WP3D = WP3DStandard;

%Following loop creates a list of latitudes and longtitudes for each
%timestep
WayPoint = [];
for k = 1:n+1
    pos = gcwaypts(pointIdentifiers(k,1),pointIdentifiers(k,2),...
        pointIdentifiers(k+1,1),pointIdentifiers(k+1,2),...
        ceil(pointIdentifiers(k+1,4)/TimeResolution)-...
        ceil(pointIdentifiers(k,4)/TimeResolution));
    WayPoint = [WayPoint, [pos(2:end, 1)';pos(2:end, 2)']];
end

%Adds 360 degress to the positions to bring to a continuous coordinate
%system for ease of computation
indices = find(WayPoint(2, :) < 0);
WayPoint(2, indices) = WayPoint(2, indices) + 360;

%Following loop creates a listing og course over ground and modifies the
%WP3D matrix so the sampled route in time space is created.
COG = zeros(1, length(WayPoint(1, :)));
for j = 1:length(WayPoint(1, :)) - 1
    dx = WayPoint(1, j+1) - WayPoint(1, j);
    dy = WayPoint(2, j+1) - WayPoint(2, j);
    if dx == 0
        if dy == 0
        elseif dy > 0
            COG(j) = 90;
        else
            COG(j) = 270;
        end
    elseif dx > 0
        if dy >= 0
            COG(j) = atand(dy/dx);
        else
            COG(j) = atand(dy/dx) + 360;
        end
    elseif dx < 0
        if dy == 0
            COG(j) = 180;
        else
            COG(j) = atand(dy/dx) + 180;
        end
    end
    
    latitudes = (linspace(WayPoint(1, j), WayPoint(1, j+1), npoints));
    longitudes = (linspace(WayPoint(2, j), WayPoint(2, j+1), npoints));
    for h=1:npoints
        latnow = round(((round(latitudes(h)*Resolution)/Resolution)-...
            f1(1))*Resolution);
        lonnow = round(((round(longitudes(h)*Resolution)/Resolution)-...
            f2(1))*Resolution);
        WP3D(latnow, lonnow, j) = 1;
    end
end
